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Three fermions with strongly repulsive interactions in a spherical harmonic trap, constitute the 
simplest nontrivial system that can exhibit the onset of itinerant ferromagnetism. Here, we present 
exact solutions for three trapped, attractively interacting fermions near a Feshbach resonance. We 
analyze energy levels on the upper branch of the resonance where the atomic interaction is effec- 
tively repulsive. When the s-wave scattering length a is sufflciently positive, three fully polarized 
fermions are energetically stable against a single spin-flip, indicating the possibility of itinerant 
ferromagnetism, as inferred in the recent experiment. We also investigate the high-temperature 
thermodynamics of a strongly repulsive or attractive Fermi gas using a quantum virial expansion. 
The second and third virial coefficients are calculated. The resulting equations of state can be tested 
in future quantitative experimental measurements at high temperatures and can provide a useful 
benchmark for quantum Monte Carlo simulations. 



I. INTRODUCTION 

Few-particle systems have become increasingly crucial 
to the physics of strongly interacting ultracold quantum 
gases [I'-'S]. Because of large interaction parameters, con- 
ventional perturbation theory approaches to quantum 
gases such as mean-field theory simply break down 0tJ| ■ 
A small ensemble of a few fermions and/or bosons, which 
is either exactly solvable or numerically tractable, is more 
amenable to nonperturhative quantal calculations. Al- 
though challenging experimentally, such ensembles bene- 
fit from the same unprecedented controllability and tun- 
ability as in a mesoscopic system containing a hundred 
thousand particles. The atomic species, the quantum 
statistics, the s-wave and higher partial wave interac- 
tions [5|, and the external trapping environment can all 
be controlled experimentally. The study of few-particle 
systems can therefore give valuable insights into the more 
complicated mesoscopic many-body physics of a strongly 
interacting quantum gas. In addition to qualitative in- 
sights, these solutions have already proved invaluable in 
developing high-temperature quantum virial or cluster 
expansions for larger systems Ja|, which have been re- 
cently verified experimentally [7|. 

The purpose of this paper is to add a further milestone 
in this direction. By exactly solving the eigenfunctions 
of three attractively interacting fermions in a spherical 
harmonic trap, we aim to give a few-body perspective of 
itinerant ferromagnetism in an effectively repulsive Fermi 
gas, which was observed as a transient phenomenon in a 
recent measurement at MIT [8j. This is possible because 
the quantum three-body problem with s-wave interac- 



* Electronic address: I xiajiliu@swin.e du. au| 
^Electronic address: hhu@swin.edu. au] 
■••Electronic address: pdrumniond@swin.edu. au| 



tions is exactly soluble in three dimensions. It is inter- 
esting to recall that the corresponding classical three- 
body problem is notoriously insoluble. The reason for 
this unexpectedly docile quantum behaviour is that the 
s-wave interaction Hamiltonian applicable to ultra-cold 
Bose and Fermi gases is essentially just a boundary condi- 
tion on an otherwise non-interacting quantum gas. Thus, 
we have an unusual situation where quantum mechanics 
actually simplifies an intractable classical problem. 

For Bose-Einstein condensates (BEC) in the strongly 
interacting regime, three trapped bosonic atoms with 
large s-wave scattering length were already investigated 
theoretically as a minimum prototype [9] of this few- 
body physics. To understand the fascinating crossover 
from a BEC to a Bardeen-Cooper-Schrieffer (BCS) super- 
fluid, two spin-up and two spin-down fermions in a trap 
were also simulated numerically, constituting the sim- 
plest model of the BEC-BCS crossover problem [1Q|, [11| . 
Moreover, knowledge of few-particle processes such as 
three-body recombination is primarily responsible for 
controlling the loss rate or lifetime of ultracold atomic 
gases, which, in many cases, imposes severe limitations 
on experiments. Important examples in this context in- 
clude the confirmation of stability of dimers in the BEC- 
BCS crossover [12| and the discovery of the celebrated 
Efimov state (i.e., a bound state of three resonantly inter- 
acting bosons) as well as the related universal four-body 
bound state (l3| . 

Whether an itinerant Fermi gas with repulsive inter- 
actions exhibits ferromagnetism is a long-standing prob- 
lem in condensed matter physics [1J|. It has recently at- 
tracted increasing attention in the cold-atom community 
[15-24] . The answer depends on a competition between 
the repulsive interaction energy and the cost of kinetic 
energy arising from Pauli exclusion. A strong repulsive 
interaction can induce polarization or ferromagnetism, 
since fermions with the same spin orientation are pro- 
tected from local interactions by the exclusion principle. 




Figure 1: (Color online) Energy spectrum of the relative mo- 
tion of a trapped two-fermion system near a Fesiibach res- 
onance (i.e, dja — 0, where d is the characteristic harmonic 
oscillator length). For a positive scattering length a > in the 
right part of the figure, the ground state is a molecule with size 
a, whose energy diverges as Erei — —ti^/{ma^). The excited 
states or the upper branch of the resonance may be viewed 
as the Hilbert space of a "repulsive" Fermi gas with the same 
scattering length a. In this two-body picture, the level from 
the point 2 to 3 is the ground state energy level of the repul- 
sive two-fermion sub-space, whose energy initially increases 
linearly with increasing a from 1.5/kj at the point 2 and fi- 
nally saturates towards 2.5hu} at the resonance point 3. For 
comparison, we illustrate as well the ground state energy level 
in the case of a negative scattering length and show how the 
energy increases with increased scattering length from point 
1 to 2. 



This, however, increases the Fermi energy, as all fermions 
must now occupy the same band. The difficulty in finding 
the transition point is that quantum correlations change 
the interaction energy in a way that is difficult to calcu- 
late in general. 

As early as the 1930's, Stoner showed with a sim- 
ple model using mean-field theory that ferromagnetism 
in a homogeneous Fermi gas will always take place 
[lj|. This model, however, gives the unphysical result 
that the interaction energy within the mean-field ap- 
proximation scales linearly with the s-wave scattering 
length a and therefore could be infinitely large. The 
predicted critical interaction strength at zero tempera- 
ture, {kpa)c = 7r/2, where kp is the Fermi wave-vector 
is also too large in the mean-field picture. An im- 
proved prediction from second-order perturbation the- 
ory, (kFa)c — 1.054, suffers from similar doubts about 
its validity [1^, [la|. Most recently, three independent 
ab-initio quantum Monte Carlo simulations conclusively 
reported a zero-temperature ferro mag netic phase transi- 
tion at ikFa)c ^ 0.8 - 0.9 [H, [H, [231 1. 

Several important issues are still open, including the 
nature of transition at finite temperatures. The unitarity 
limited interaction energy at infinitely large scattering 
length (a — > c») is also to be determined. 



The exciting experiment at MIT of ^Li atoms realized 
in some sense the textbook Stoner model Q. A crucial 
aspect in this realization is that the interatomic interac- 
tions are very different from the conventional model of 
hard-sphere interactions 1 15l. \22L |23| . In the experiment, 
the atoms are on the upper branch of a Feshbach reso- 
nance with a positive scattering length a > and a neg- 
ligible effective interaction range. The properties of the 
atoms are therefore universal, independent of the details 
of the interactions [2^ [2g|. This universality, however, 
comes at a price: the underlying two-body interaction 
is always attractive, so that the ground state is a gas of 
molecules of size a. The experiment thus suffers from 
considerable atom loss and has to be carried out under 
nonequilibirum conditions. This is clearly explained in 
Fig. 1 , which shows the relative energy spectrum of a pair 
of fermions in a harmonic trap with frequency uj across a 
Feshbach resonance [27]. For a > 0, the whole spectrum 
consists of two distinct parts, the lowest ground-state 
branch diverging as Erei — —h^/{ma'^) and the regular 
upper branch having a finite energy. In the context of 
this two-body picture, an s-wave "repulsive" Fermi gas 
is realized, provided that there are no pairs of fermions 
occupying the ground state branch of molecules. How- 
ever, as far as the many-particle aspect is concerned, it 
is not clear to what extent this two-body picture of a 
"repulsive" Fermi gas will persist. In other words, can we 
prove rigorously that tlie whole Hilbert space of an at- 
tractive many-fermion system with a positive scattering 
length consists exactly of a sub- Hilbert space of a repul- 
sive Fermi gas with the same scattering length, together 
with an orthogonalized subspace of molecules? 

This paper addresses the problem of itinerant ferro- 
magnetism in an attractive Fermi gas using a few-particle 
perspective, by examining the exact solutions for the en- 
ergy spectrum of three trapped, attractively interacting 
fermions in their upper branch of the Feshbach resonance. 
Our main results may be summarized as follows. 

• First, we present an elegant and physically trans- 
parent way to exactly solve the Hamiltonian of 
three interacting fermions in a harmonic trap. The 
method may easily be generalized to treat other 
systems with different types of atomic species, ge- 
ometries, and interactions. 

• Secondly, we observe clearly from the whole relative 
energy spectrum of three attractive fermions (see 
Fig. 3) that, there are indeed two branches of the 
spectrum on the side of positive scattering length. 
As the scattering length goes to an infinitely small 
positive value, the lower branch diverges in energy 
to — oo, while another upper branch always con- 
verges to the non-interacting limit. The latter may 
be interpreted as the energy spectrum of three "re- 
pulsively" interacting fermions. However, close to 
the Feshbach resonance, there are many nontriv- 
ial avoided crossings between two types of spec- 
trum, making it difficult to unambiguously iden- 



tify a repulsive Fermi system. These avoided cross- 
ings are expected in more general cases, and lead to 
nontrivial consequences in a time-dependent field- 
sweep experiment passing from the weakly inter- 
acting regime at a = 0+ to the unitarity limit at 
a = +00. 

• Thirdly, we show exactly that near the Feshbach 
resonance, three "repulsively" interacting fermions 
in a spherical harinonic trap, say, in a two spin-up 
and one spin-down configuration, are higher in to- 
tal energy than three fully polarized fermions (see 
the ground state energy in Fig. 3b). Thus, there 
must be a ferromagnetic transition occurring at a 
certain critical scattering length. Note that fer- 
romagnetism cannot be obtained in a two-fermion 
system. As shown in Fig. 1, even at resonance 
the total ground state energy of a repulsively in- 
teracting pair, Epair — 4:hLu, which is the sum 
of the relative ground state energy Erei — 2.bfvjj 
and the zero-point energy of center-of-mass mo- 
tion Ecm = l-5Huj, cannot be larger than the total 
ground state energy of two fully polarized fermions, 



i.e., E^ 



n 



ATujj. 



• Last but most importantly, we obtain the high tem- 
perature equations of state of strongly interacting 
Fermi gases (see Figs. 4, 5, 6, and 7), within a quan- 
tum virial expansion theory, which was developed 
recently by the present authors [6] . The second and 
third virial (expansion) coefficients of both attrac- 
tive and repulsive Fermi gases can be calculated, 
using the full energy spectrum of three interacting 
fermions. In the unitarity limit, we find that in the 
high temperature regime where our quantum virial 
expansion is reliable, the itinerant ferromagnetism 
disappears. 

The paper is organized as follows. In the next section, 
we outline the theoretical model for a few fermions with 
s-wave interactions in a spherical harmonic trap. In Sec. 
Ill, we explain how to construct the exact wavefunctions 
for three interacting fermions and discuss in detail the 
whole energy spectrum. In Sec. IV we develop a quan- 
tum virial expansion for thermodynamics and calculate 
the second and third virial coefficients, based on the full 
energy spectrum of two-fermion and three-fermion sys- 
tems, respectively. The high temperature equations of 
state of strongly interacting Fermi gases are then calcu- 
lated and discussed in Sec. V. Finally, Sec. VI is devoted 
to conclusions and some final remarks. The appendix 
shows the numerical details of the exact solutions. 



II. MODELS 

Consider a few fermions in a three-dimensional (3D) 
isotropic harmonic trap V{'x) = muP'x^ jl with the same 
mass m and trapping frequency w, occupying two dif- 
ferent hyperfine states or two spin states. The fermions 



with unlike spins attract each other via a short-range 
s-wave contact interaction. It is convenient to use the 
Bethe-Peierls boundary condition to replace the s-wave 
pseudopotential. That is, when any particles i and j with 
unlike spins close to each other, rij — |x,; — x^j — >■ 0, the 
few-particle wave functio n •i/ > (xi ,X2, ...,xjv) with proper 
symmetry should satisfy [28|-|30|, 



i) ^ Aij (Xij = 






(1) 



ij, 



where Aij(X.ij, {^k^ij}) is a function independent of . 
and a is s-wave scattering length. This boundary condi- 
tion can be equivalently written as. 



lim 



d (rjj-0) _ njTp 
drji a 



(2) 



Otherwise, the wave function ip obeys a non-interacting 
Schrodinger equation. 
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We now describe how to solve all the wave functions with 
energy level E for a two- or three-fermion system. 



III. METHOD 

In a harmonic trap, it is useful to separate the center- 
of-mass motion and relative motion. We thus define the 
following center-of-mass coordinate R and relative coor- 
dinates Vi (i > 2) for N fermions in a harmonic trap 

MM, 



R 



xi 



xjv) /N, 



and 



r, = 



i- 1 



(""^H 



(4) 



(5) 



respectively. In this Jacobi coordinate, the Hamiltonian 
of the non-interacting Schrodinger equation takes the 
form Hq = Han + 'Hrei, whcrc. 
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The center-of-mass motion is simply that of a harmoni- 
cally trapped particle of mass A/ = Nm, with well-known 
wave functions and spectrum Ecm — [ncm + S/2)huj, 
where ricm = 0,1,2... is a non-negative integer. In the 
presence of interactions, the relative Hamiltonian should 
be solved in conjunction with the Bethe-Peierls boundary 
condition, Eq. ^. 



A. Two fermions in a 3D harmonic trap 

Let us first briefly revisit tlie two-fermion problem in a 
harmonic trap, where the relative Schrodinger equation 
becomes 



-2^1^^ 
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V2t{v)^ErelV2t(j), 



(8) 



where two fermions with unlike spins do not stay at the 
same position (r > 0). Here, we have re-defined r = \/2t2 
and have used a reduced mass /i — m/2. It is clear that 
only the I = subspace of the relative wave function 
is affected by the s-wave contact interaction. According 
to the Bethe-Peierls boundary condition, as r — ?► the 
relative wave function should take the form, 'il^2f{r) — ?> 
(1/r- 1/a), or satisfy, d {rip^2b) jdr = - {r^lf) I a. The 
two-fermion problem in a harmonic trap was first solved 
by Busch and coworkers [27]. Here, we present a simple 
physical interpretation of the solution. 

The key point is that, regardless of the boundary con- 
dition, there are two types of general solutions of the 
relative Schrodinger equation ([5]) in the I = subspace, 
^2b'(7') ex exp(— r^/2d^)/(r/d). Here the function ]{x) 
can either be the first kind of Kummer confluent hyper- 
geometric function \F\ or the second kind of Kummer 
confluent hypergeometric function JJ . We have taken 
d = \Jhl [110 as the characteristic length scale of the trap. 
In the absence of interactions, the first Kummer func- 
tion gives rise to the standard wave function of 3D har- 
monic oscillators. With interactions, however, we have 
to choose the second Kummer function C/, since it di- 
verges as 1/r at origin and thus satisfies the Bethe-Peierls 
boundary condition. 

Therefore, the (un-normalized) relative wave function 
and relative energy should be rewritten as. 



i^\t{T-.v) = n-v)u{^ 
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Figure 2: (Color online) Configuration of three interacting 
fermions, two spin-up and one spin-down. 



{p-rei = 0, 1, 2...), which recovers the spectrum in the non- 
interacting limit. With increasingly attractive interac- 
tions, the energies decrease. In the unitarity (resonance) 
limit where the scattering length diverges, a — >■ ±cxi, we 
find that v(a = ±(xi) = rirei — 1/2. As the attraction 
increases further, the scattering length becomes positive 
and decreases in magnitude. We then observe two dis- 
tinct types of behavior: the ground state is a molecule of 
size a, whose energy diverges asymptotically as —h^/ma^ 
as a — >■ 0^, while the excited states may be viewed as two 
repulsively interacting fermions with the same scattering 
length a. Their energies decrease to the non- interacting 
values as a -^ 0+. 

In this two-body picture, a universal repulsively inter- 
acting Fermi gas with zero-range interaction potentials 
may be realized on the positive scattering length side of 
a Feshbach resonance for an attractive interaction po- 
tential, provided that all two fermions with unlike spins 
occupy the exited states or the upper branch of the two- 
body energy spectrum. 



and 



Erel = (2l^ + -)flhJ, 



(10) 



respectively. Here, F is the Gamma function, the 
real number u plays the role of a quantum number 
and should be determined by the boundary condition, 
Vmvr^o d {r'4>2f^ Idr = — (?'V'2b') Z*^- By examining the 
short range behavior of the second Kummer function 
C/(— i^, 3/2, x), this leads to the familiar equation for en- 
ergy levels, 



2F(-i/) _ d 
Y{-v-\l2) ~~a 



(11) 



In Fig. 1, we give the resulting energy spectrum as a 
function of the dimensionless interaction strength dj a. 

The spectrum is easy to understand. At infinitely 
small scattering length a — >■ 0~, v(a = 0^) — rirei 



B. Three fermions in a 3D harmonic trap 

Let us turn to the three fermion case by considering 
two spin- up fermions and one spin-down fermion, i.e., the 
tit configuration shown in Fig. 2. The relative Hamil- 
tonian can be written as 1291 ISOll , 



H.eZ = |^(V? + V^) + ^M-M-'+p')- (12) 

where we have redefined the Jacobi coordinates r — -\/2r2 
and p — v2r3, which measure the distance between the 
particle 1 and 2 (i.e., pair), and the distance from the 
particle 3 to the center-of-mass of the pair, respectively. 



1. General exact solutions 



Thus, the Bethe-Peierls boundary condition becomes, 



Inspired by the two-fermion solution, it is readily seen 
that the relative wave function of the Hamiltonian P^ 
may be expanded into products of two Kummer conflu- 
ent hypergeometric functions. Intuitively, we may write 
down the following ansatz [6|, 



V'3t(r,p) = (l-7'i3)x(r,p), 



(13) 



where, 



X (r, p)^Yl «"^2b' (^; ^'.")^ni (p) vr ip) . (14) 



The two-body relative wave function ipl^f (r; vi^n) with en- 
ergy {2i/i^n + i/2)huj describes the motion of the paired 
particles 1 and 2, and the wave function Rni (p) Yi™' (p) 
with energy (2n + / + 3/2)fiaj gives the motion of particle 
3 relative to the pair. Here. Rni (p) is the standard radial 
wave function of a 3D harmonic oscillator and Y"^ (p) is 
the spherical harmonic. Owing to the rotational symme- 
try of the relative Hamiltonian ((T^), it is easy to see that 
the relative angular momenta I and m are good quantum 
numbers. The value of ;//,„ is uniquely determined from 
energy conservation, 

Erei = [{2j^i.n + 3/2) + (2n + I + 3/2)] huj, (15) 

for a given relative energy Erei- It varies with the index 
n at a given angular momentum I. Finally, T'la is an 
exchange operator for particles 1 and 3, which ensures 
the correct exchange symmetry of the relative wave func- 
tion due to Fermi exclusion principle, i.e., T^isx (r, p) = 
X (r/2 + V3p/2, V3r/2 - p/2). The relative energy Erei 
together with the expansion coefficient a„ should be de- 
termined by the Bethe-Peierls boundary condition, i.e., 
\iu].r^o[dri{}lf {r,p)]/dr = — [ri/^g^' (r,p)]/a. We note 
that the second Bethe-Peierls boundary condition in case 
of particle 2 approaching particle 3 is satisfied automati- 
cally due to the exchange operator acting on the relative 
wave function. 

By writing ^(r,/?) — (p{r, p)Yf^ (p), the Bethe-Peierls 
boundary condition takes the form (r — > 0), 
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Using the asymptotic behavior of the second kind of 
Kummer function, YaUx^^^ {—vi^n)U{—i'i^n,'i/2,x'^) — 
^JH jx — 2^Y {—lyi^n) /r {—I'l.n — 1/2), it is easy to show 
that in the limit of r — >■ 0, 

--[r0(r,p)]=-^^a„i?„,(p), (17) 
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Projecting onto the orthogonal and complete set of basis 
functions Rni (p), we find that a secular equation. 
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where we have defined the matrix coefficient 
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which arises from the exchange effect due to the operator 
7^13 . In the absence of C„„' , the above secular equation 
describes a three-fermion problem of a pair and a sin- 
gle particle, un-correlated to each other. It then simply 
reduces to Eq. (|lip . as expected. 

The secular equation (PT|) was first obtained by Kestner 
and Duan by solving the three-particle scattering prob- 
lem using Green function [3l|. To solve it, for a given 
scattering length we may try different values of relative 
energy Erei, implicit via vi^n- However, it turns out to be 
more convenient to diagonalize the matrix A = {A„„'} 
for a given relative energy, where 



2r(-.,„) (-1)' 

"" "r(-^,„-i/2)'^""'+ 0F 
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(23) 



The eigenvalues of the matrix A then gives all the pos- 
sible values of d/a for a particular relative energy. We 
finally invert a{Erei) to obtain the relative energy as a 
function of the scattering length. Numerically, we find 
that the matrix A is symmetric and thus the standard 
diagonalization algorithm can be used. We outline the 
details of the numerical calculation of Eq. (P5)) in the 
Appendix A. 



2. Exact solutions in the unitarity limit 

In the unitarity limit with infinitely large scatter- 
ing length, a — > oo, we may obtain more physical so- 
lutions using hyperspherical coordinates, as shown by 
Werner and Castin |28l . [30| . By defining a hyperradius 
R = \J{r'^ + P^)/2 and hyperangles il = (a,r^p\ where 
a = arctan(r/p) and f and p are respectively the unit 
vector along r and p, we may write |28l . |30|, 



V'i 
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to decouple the motion in the hyperradius and hyperan- 
gles for given relative angular momenta / and m . It leads 
to the following decoupled Schrodinger equations [30| . 



^ + LJ^R^ F = 2ErelF, (25) 




and 






.^(«), (26) 



where sf ^ is the eigenvalue for the n-ih wave function of 
the hyperangle equation. 

For three- fermions, sf^ is always positive. There- 
fore, the hyperradius equation ([25]) can be interpreted 
as a Schrodinger equation for a fictitious particle of mass 
unity moving in two dimensions in an effective potential 
{sf^/R^ + to'^R^) with a bounded wave function F{R). 
The resulting spectrum is |28l . [3(J | 



Erei = (2g-t- s;,„ + l)hu, 



(27) 



where the good quantum number q labels the number of 
nodes in the hyperradius wave function. 

The eigenvalue Si_„ should be determined by the Bethe- 
Peierls boundary condition, which in hyperspherical co- 
ordinates takes the from [2a, [Sfl , 



(28) 
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In addition, we need to impose the boundary condition 
(f (7r/2) = 0, since the relative wave function ([M)) should 
not be singular at a = 7r/2. The general solution of the 
hyperangle equation (^5)) satisfying (p {n/2) = is given 
by, 



ifi ex x^+\Fi 



l + l-si^„ l + l 



si, 



^■' + 1 



(29) 

where x = cos(a) and 2F1 is the hypergeometric function. 
In the absence of interactions, the Bethe-Peierls bound- 
ary condition ([28| should be replaced by ip (0) = 0, since 
the relative wave function (|24p should not be singular 
at a = either. As ^(0) = r{l + 3/2)r(l/2)/[r((/ + 
2 + •Si,n)/2)r((Z + 2 — sin)/2)], this boundary condition 
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leads to [^ + 2 

where n — 0, 1, 2, ... is a non-negative integer and we have 
used the superscript "1" to denote a non- interacting sys- 
tem. However, a spurious solution occurs when I = and 
n = 0, for which sj „ = 2, 93(a) — sin(2a)/2 and thus, the 
symmetry operator (1 — V13) gives a vanishing relative 
wave function in Eq. (PH) that should be discarded |30| . 
We conclude that for three non-interacting fermions. 



,(1) 



2n + 4, 1 = 
2n + l + 2, l>0 



(30) 



For three interacting fermions, we need to determine si^n 
by substituting the general solution ([^^ into the Bethe- 
Peierls boundary condition (pS)) . In the Appendix B, we 
describe how to accurately calculate si^n- In the bound- 
ary condition Eq. (P5)l . the leading effect of interactions 
is carried by ip' (0) and therefore, ip' (0) = determines 
the asymptotic values of s;^„ at large momentum I or n. 
This gives rise to (/ -f 1 — Si_„)/2 = —n, or. 



Sin 



2n + 3, 1 = 
2n + l + l, l> 



(31) 



where we have used a bar to indicate the asymptotic re- 
sults. By comparing Eqs. ([50)) and ([5T|) . asymptotically 
the attractive interaction will reduce s;_„ by a unity. 



3. 



Energy spectrum of three interacting fermions 



We have numerically solved both the general exact so- 
lution ([T5| along the BEC-BCS crossover and the exact 
solution (j24p in the unitarity limit. In the latter uni- 
tary case, the accuracy of results can be improved to 
arbitrary precision by using suitable mathematical soft- 
ware, described in Appendix B. Fig. 3 reports the energy 
spectrum of three interacting fermions with increasingly 
attractive interaction strength at three total relative an- 
gular momenta, Z = 0, 1, and 2. For a given scattering 
length, we typically calculate several ten thousand en- 
ergy levels (i.e., Erei < {I + 256) fiw) in each different 
subspace. To construct the matrix A, Eq. ([^5]) . we have 
kept a maximum value of rimax — 128 in the functions 
Rni (p). Using the accurate spectrum in the unitarity 
limit as a benchmark, we estimate that the typical rela- 
tive numerical error of energy levels is less than 10~^. We 
have found a number of nontrivial features in the energy 
spectrum. 

The spectrum on the BCS side is relatively simple. 
It can be understood as a non-interacting spectrum at 
d/a -^ —00, in which Erei — {2Q + 3)huj at I = and 
Erei — (2Q + I + l)huj at / > 1, with a positive integer 
Q = 1,2,3,... that denotes also the degeneracy of the 
energy levels. The attractive interactions reduce the en- 
ergies and at the same time lift the degeneracy. Above 
the resonance or unitary point of d/a = 0, however, the 
spectrum becomes much more complicated. 

There are a group of nearly vertical energy levels that 
diverge towards the BEC limit of d/a — >■ -l-oo. From the 
two-body relative energy spectrum in Fig. 1, we may 
identify these as energy states containing a molecule of 
size a and a fermion. For a given scattering length, these 
nearly vertical energy level differ by about 2huj, resulting 
from the motion of the fermion relative to the molecule. 
In addition to the nearly vertical energy levels, most in- 
terestingly, we observe also some nearly horizontal energy 
levels, which converge to the non-interacting spectrum in 
the BEC limit. In analogy with the two-body case, we 
may identify these horizontal levels as the energy spec- 
trum of three repulsively interacting fermions. We show 




Figure 3: (Color online) Relative energy spectrum of three 
interacting fermions at different subspaces or relative angular 
momenta I. On the positive scattering length (BEC) side of 
the resonance, there are two types of energy levels: one (is 
vertical and) diverges with decreasing the scattering length a 
and the other (is horizontal) converges to the non-interacting 
spectrum. The latter may be viewed as the energy spectrum 
of three repulsively interacting fermions. In analogy with the 
two-fermion case, we show in the ground state subspace {I — 
1) the ground state energy level of the repulsive three-fermion 
system (i.e, from point 2 to 3), as well as the ground state 
energy level of the attractive three-fermion system for a < 
(i.e., from the point 1 to 2). In the unitarity limit, we show 
by the circles the energy levels that should be excluded when 
we identify the energy spectrum for infinitely large repulsive 
interactions. 



explicitly in Fig. 3b the ground state level of three re- 
pulsively interacting fermions, which increases in energy 
from the point 2 to 3 with increasing scattering length 
from a = 0"*" to a = -l-oo. For comparison, we also show 
the ground state level of three attractively interacting 
fermions at a negative scattering length, which decreases 
in energy from the point 2 to 1 with increasing absolute 
value of a. 

This identification of energy spectrum for repulsive in- 
teractions, however, is not as rigorous as in the two-body 
case. There are many apparent avoided crossings be- 
tween the vertical and horizontal energy levels. There- 
fore, by changing a positive scattering length from the 
BEC limit to the unitarity limit, three fermions initially 
at the horizontal level may finally transition into a verti- 
cal level, provided that the sweep of scattering length is 
sufficiently slow and adiabatic. This leads to the conver- 
sion of fermionic atoms to bosonic molecules. A detailed 
analysis of the loss rate of fermionic atoms as a func- 
tion of sweep rate may be straightforward obtained by 
applying the Landau-Zener tunnelling model. 

Let us now focus on the resonance case of most signifi- 
cant interest. In Fig. 3, we show explicitly by green dots 
the vertical energy levels in the unitarity limit. These 
levels should be excluded if we are interested in the spec- 
trum of repulsively interacting fermions. Amazingly, for 
each given angular momentum, these energy levels form 
a regular ladder with an exact energy spacing 2fiuj [29|. 
Using the exact solution in the unitarity limit, Eq. (|27p . 
we may identify unambiguously that the energy ladder is 
given by. 



Erel = {2q + S;,o -|- 1) fuxJ . 



(32) 



Therefore, in the unitarity limit the lowest-order solution 
of the hyperangle equation gives rise to the relative wave 
function of a molecule and a fermion. Thus, it should be 
discarded when considering three resonantly interacting 
fermions with an effective repulsive interaction. 

This observation immediately leads to the ground state 
energy of three repulsively interacting fermions, 

£;tit ^ (s^_^ + l)huj + l.5huj ~ 6.858249309fta;, (33) 

including the zero-point energy of the center-of-mass mo- 
tion, 1.5haj. This ground state energy is higher than that 
of three fully polarized fermions, which is. 



4r = i-5/^- 



2.5huj + 2.5hu} = 6.5hu}. 



(34) 



Thus, in the presence of repulsive interactions, the 
ground state of three fully polarized fermions is unsta- 
ble with respect to a single spin-fiip, suggestive of an 
itinerant ferromagnetic transition at a certain scattering 
length for three fermions. 



IV. QUANTUM VIRIAL EXPANSION FOR 
THERMODYNAMICS 



The few-particle solutions presented above can provide 
information about the high temperature thermodynamics 
of many-body systems, through a quantum virial expan- 
sion of the grand thermodynamic potential [6, 32'| . In the 
grand canonical ensemble, the thermodynamic potential 
is given by, 

n^-ksTlnZ, (35) 

where kB is the Boltzmann constant and 

Z = Tr exp [- (H - /lA/") /ksT] (36) 

is the grand partition function. We may rewrite this in 
terms of the partition function of clusters. 



Q„ = Tr„ [exp (-H/fcsT)] 



(37) 



where the integer n denotes the number of particles in the 
cluster and the trace Tr„ is taken over n-particle states 
with a proper symmetry. The partition function of clus- 
ters Qn can be calculated using the complete solutions 
of a n-particles system. The grand partition function is 
then written as 



Z = l + zQi+z'^Q2 + 



(38) 



where z = exp {^./ksT) is the fugacity. At high tem- 
peratures, it is well-known that the chemical potential 
H diverges to — oo, so the fugacity would be very small, 
z <C 1. We can then expand the high-temperature ther- 
modynamic potential in powers of the small parameter 



A. Non-interacting virial coefficients 

The background non-interacting virial coefficients 
can be conveniently determined by the non-interacting 
thermodynamic potential. For a homogeneous two- 
component Fermi gas, this takes the form, 

oo 

^11 = -^^ ;! / ^"' '- (1 + --') d^^ (44) 



where A = \^'K'h? / (jnkBT)]^/'^ is the thermal wavelength 
and Qi.hom = 2P^/A'^. This leads to 



,(1) _ (-1) 



n+1 



.horn. 



,5/2 • 



(45) 



Hereafter, we use the subscript "hom" to denote the quan- 
tity in the homogeneous case, otherwise, by default we 
refer to a trapped system. 

For a harmonically trapped Fermi gas, the non- 
interacting thermodynamic potential in the semiclassical 
limit (neglecting the discreteness of the energy spectrum) 
is. 



{hwf 2 



t^\a{l + ze'-^)dt, (46) 



where Qi = 2 (ksT) / [hw) . Taylor-expanding in pow- 
ers of z gives rise to 



hi 



1) _ ("1) 



n+l 



(47) 



Q. = -kBTQi [z + b2z^ + ■■■ + 6„z" + •••], (39) 

where 6„ may be referred to as the n-th (virial) expansion 
coefficient. It is readily seen that, 

&2 - {Q2-Ql/2)/Qi, (40) 

bs = (g3-QiQ2 + Q?/3)/Qi, etc. (41) 

These equations present a general definition of the quan- 
tum virial expansion and are applicable to both homo- 
geneous and trapped systems. The determination of the 
n-th virial coefficient requires knowledge of up to the n- 
body problem. 

In practice, it is convenient to concentrate on the in- 
teraction effects only. We thus consider the difference 
A6„ = 6„ — 6„ and AQn = Qn — Qn , where the super- 
script "1" denotes the non-interacting systems. For the 
second and third virial coefficient, we shall calculate 



and 



A&2 = AQ2/Q1 



A63 = Ag3/Qi-AQ2. 



(42) 



(43) 



We note that the non-interacting virial coefficients in the 
homogeneous case and trapped case are related by. 



7 (l^ ?T,,hoin 



(48) 



B. Second virial coefficient in a harmonic trap 

We now calculate the second virial coefficient of a 
trapped interacting Fermi gas. In a harmonic trap, the 
oscillator length d provides a large length scale, compared 
to the thermal wavelength A. Alternatively, we may use 
Lj = huj/kBT <C 1 to characterize the intrinsic length 
scale relative to the trap. All the virial coefficients and 
cluster partition functions in harmonic traps therefore 
depend on the small parameter iJj. We shall be interested 
in a universal regime with vanishing w, in accord with 
the large number of atoms in a real experiment. 

To obtain A&2, we consider separately AQ2 and Qi. 
The single-particle partition function Qi is determined 
by the single-particle spectrum of a 3D harmonic os- 
cillator, Erd = (2n + l + 3/2)hu}. We find that Qi = 
2/[exp(+w/2) - exp(-w/2)]3 ~ 2{kBTf / {fuvf. The 



prefactor of two accounts for the two possible spin states 
of a single fermion. In the calculation of IS.Q2 , it is easy 
to see that the summation over the center-of-mass energy 
gives exactly Qi/2. Using Eq. HU]), we find that, 



A5, * = 



\n 



-(2iy„+3/2)(i _ -(2iy^^'+3/2)ii 



(49) 



To proceed, it is important to analyze analytically the 
behavior of iJ^eZ at high energies. For this purpose, we 
introduce a relative energy E^-ei , which is the solution of 
Eq. (P5)l in the absence of the exchange term Cmn, and 
can be constructed directly from the two-body relative 
energy. In the subspace with a total relative momentum 
L it takes the form, 



where the non-interacting i^n — n [n ~ 0,1,2,...) and 
the superscript "att" (or "rep") means the coefficient of an 
attractively (or repulsively) interacting Fermi gas. The 
second virial coefficient of a trapped attractive Fermi gas 
in the BEC-BCS crossover was given in Fig. 3a of Ref. 

To consider the second virial coefficient of a repulsively 
interacting Fermi gas, we shall restrict ourselves to a posi- 
tive scattering length and exclude the lowest ground state 
energy level in the summation of the first term in Eq. 
(|^ . which corresponds to a bound molecule. 



1. Unitarity limit 

At resonance with an infinitely large scattering length, 
the spectrum is known exactly: i'n,oo = n — 1/2, giving 
rise to. 



Auatt 
^''2.cx) 



1 exp {-Cj/2) _ 1 J_ 

2 [1 -t- exp (-w)] ^ ^4 " 32^ 



~2 



(50) 



For a repulsive Fermi gas in the unitarity limit, we shall 
discard the lowest 'molecular' state with 1/9,00 = ~l/2 
and therefore, 



. ,rep ^ 1 exp(-a;/2) 
2'°° 2 [1 -t- exp i+^)] 



4 



(51) 



The term uj'^ or w in Eqs. (j50p and (j5ip is nonuniversal 
and is negligibly small for a cloud with a large number 
of atoms. We therefore obtain the universal second virial 
coefficients: A6^*^ = -1-1/4 and Ab^^;^ = -1/4, which 
are temperature independent. 



Erei = {2n + l + 3/2) hLj + {2i^ + 3/2)fiw, 



(53) 



where v is the solution of the two-body spectrum of Eq. 
PT|) . At high energies the full spectrum Erei approaches 
asymptotically to Erei as the exchange effect becomes 
increasingly insignificant. There is an important excep- 
tion, however, occurring at zero total relative momentum 
1 = 0. As mentioned earlier, the solution of Erei at n = 
and / = is spurious and does not match any solution of 
Erei ■ Therefore, for the I = subspace, we require n > 1 
in Eq. (P)). 

It is easy to see that if we keep the spurious solution in 
the I = subspace, the difference [^exp(— i^rez/fesT") — 
^exp(— i^^^gj/fc^T)] is exactly equal to AQ2, since in 
Eq. ([55)1 the first part of spectrum is exactly identical 
to the spectrum of center-of-mass motion. The spurious 
solution gives a contribution. 



U' 



-{2i^n+3)uj 



-(2v 



= 2e-^'^/^Abf\ 



(54) 



which should be subtracted. Keeping this in mind, we 
finally arrive at the following expression for the third 
virial coefficient of a trapped Fermi gas with repulsive 
interactions: 



Abf' 



E 



Erei 



2e~^"/2A6f*. 



(55) 



The summation is over all the possible relative energy 
levels Erei and their asymptotic values Erei- It is well- 
behaved and converges at any scattering length. The 
third virial coefficient of a trapped attractive Fermi gas 
in the BEC-BCS crossover was given in Fig. 3b of Ref. 



a- 



C. Third virial coefficient in a harmonic trap 



1. Unitarity limit 



The calculation of the third virial coefficient, which is 
given by A63 = AQ3/Q1 — AQ2, is more complicated. 
Either the term AQ3/Q1 or AQ2 diverges as a; — )• 0, but 
the leading divergences cancel with each other. In the 
numerical calculation, we have to carefully separate the 
leading divergent term and calculate them analytically. 
It is readily seen that the spin states of fit and |ti 
configurations contribute equally to Q3. The term Qi in 
the denominators is canceled exactly by the summation 
over the center-of-mass energy. We thus have 



AQ3/Q1 = [Y,eM~Erel/kBT)- 



-J2eM~El'J/kBT)]. 
(52) 



In the unitarity limit, it is more convenient to use the 
exact spectrum given by Eq. (|27p . where s;_„ can be 
obtained numerically to arbitrary accuracy and the non- 
interacting s\J^ is given by Eq. ((50)) . To control the 
divergence problem, we shall use the same strategy as 
before and to approach si^n by using its asymptotic value 
s/,„ given in Eq. ^I^. 

Integrating out the q degree of freedom and using Eq. 
([50)1 to calculate AQ2; we find that. 



Auatt 



1-e- 



E 



A 



(56) 
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where A is given by 



(1-e-") 



2 • 



(57) 



We note that for the summation, imphcitly there is a 
prefactor {21 + 1), accounting for the degeneracy of each 
subspace. The value of A can then be calculated analyt- 
ically, leading to. 



A = -e-^ (1 - e-") . (58) 

= 512 and I < 



We have calculated numerically J^i „(£ 
by imposing the cut-offs of n < n,, 
^max = 512. We find that. 



\uatt 

^"3,00 



-0.06833960- 



0.038867w2 



(59) 



The numerical accuracy can be further improved by suit- 
ably enlarging rimax and Zmax- For a Fermi gas with in- 
finitely large repulsions, we need to exclude the states in- 
volving a molecule. Thus, in the calculation of AQ3/Q1, 
we exclude the energy levels associated with s;.„=Oj as 
given by Eq. (|5^ . In the calculation of AQ2, we shall 
remove the lowest two-body state with fo,oo — ~l/2- In 
the end, we find that. 



Abl% ~ 0.34976 



0.77607(1; 



(60) 



By neglecting the dependence on a) in the thermodynamic 
limit, we obtain the universal third virial coefficients: 



^"3,00 



-0.06833960, 



A6S;^ ~ 0.34976. 



(61) 
(62) 



D. Unitary virial coefHcients in homogeneous space 

We have so far studied the virial coefficients in a har- 
monic trap. In the unitarity limit, there is a simple re- 
lation between the trapped and homogeneous virial co- 
efficient, as inspired by relation (pS)) . This stems from 
the universal temperature independence of all virial co- 
efficients in the unitarity limit. In the thermodynamic 
limit, let us consider the thermodynamic potential of a 
harmonic trapped Fermi gas in the local density approx- 
imation il = J drn{r), where il(r) is the local thermo- 
dynamic potential 



rj(r) oc z (r) -I- 62,oo,hom2^ (r) + ■ 



(63) 



Here, the local fugacity z (r) — zexp[—V{r)/kBT] is de- 
termined by the local chemical potential /x(r) = ii—V{r). 
On spatial integration, it is readily seen that the univer- 
sal (temperature independent) part of the trapped virial 
coefficient is. 



'^n,c 



bn. 



cxD.hom 



l3/2 



(64) 



We therefore immediately obtain that the homogeneous 
second virial coefficients in the unitarity limit are: 





(65) 


A ,rep _ 
""2, 00, horn /K' 


(66) 


igeneous third virial coefficients are: 




^bf^Mm - -0.35501298, 


(67) 


A6S:^,ho„. - +1-8174. 


(68) 



The homogeneous virial coefficients are therefore signifi- 
cantly larger than their trapped counterparts. The factor 
of n^/^ is clearly due to the higher density of states in a 
harmonically trapped geometry. 



V. HIGH-T EQUATION OF STATE OF A 
STRONGLY INTERACTING FERMI GAS 

We are now ready to calculate the equation of states 
in the high temperature regime, by using the thermody- 
namic potential 



fih 



om — "horn 



V 



IkeT 
A3 



(A5: 



'2, ho 



.Z^ + 



and 



n = Q} 



2 [kBTY 

{hwf 



(Ab2z^ + Absz^ 



(69) 



(70) 



respectively, for a homogeneous or a harmonically 
trapped Fermi gas. Here, the non-interacting thermo- 
dynamic potentials are given by Eqs. ([^^ and (P5|. 
All the other thermodynamic quantities can be derived 
from the thermodynamic potential by the standard ther- 
modynamic relations, for example, N — —dQ/dfi, S = 
-dfl/dT, and then E = n + TS + ^iN. 

As an concrete example, let us focus on the unitarity 
limit in the thermodynamic limit, which is of the great- 
est interest. The equations of state are easy to calculate 
because of the temperature independence of virial coef- 
ficients. It is also easy to check the well-known scaling 
relation in the unitarity limit: E = —3fl/2 for a homo- 
geneous Fermi gas 1 321] and E = — 3il for a harmonically 
trapped Fermi gas Q ■ The difference of the factor of two 
arises from the fact that according to the virial theorem, 
in harmonic traps the internal energy is exactly equal to 
the trapping potential energy. 

To be dimensionless, we take the Fermi temperature 
Tp or Fermi energy (Ep — ksTp) as the units for tem- 
perature and energy. For a homogeneous or a harmon- 
ically trapped Fermi gas, the Fermi energy is given by 
Ef = ft2(37r2jV/T/)2/3/2TO and Ep = (37V)i/^fia;, respec- 
tively. In the actual calculations, we determine the num- 
ber of atoms N , the total entropy S, and the total energy 
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E at given fugacity and a fixed temperature (i.e., T — 1), 
and consequently obtain the Fermi temperature Tp and 
Fermi energy Ep ■ We then plot the energy or energy per 
particle, E/{NEf) and S/{NkB), as a function of the 
reduced temperature T /Tp. 




Figure 4: (Color online) Energy per particle E/{NEf) as 
a function of reduce temperature T/Tf for a homogeneous 
Fermi gas with infinitely attractive and repulsive interactions. 
The predictions of quantum virial expansion up to the second- 
and third-order are shown by solid line and dashed line, re- 
spectively. For comparison, we plot the ideal gas result by 
the dot-dashed line. 



A. Homogeneous equation of state 

We report in Figs. 4 and 5 the temperature depen- 
dence of energy and entropy of a strongly attractively or 
repulsively interacting homogeneous Fermi gas. The solid 
line and dashed line are the predictions of the quantum 
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Figure 5: (Color online) Entropy per particle S/{NEf) as 
a function of reduce temperature T/Tf for a homogeneous 
Fermi gas with infinitely attractive and repulsive interactions. 
The others are the same as in Fig. 4. 




T/T 



Figure 6: (Color online) Energy per particle E/{NEf) as a 
function of reduced temperature T/Tf for a trapped Fermi 
gas with infinitely attractive and repulsive interactions. The 
predictions of quantum virial expansion up to the second- and 
third-order are shown by solid line and dashed line, respec- 
tively. For comparison, we plot the ideal gas result by the 
dot-dashed line. We show also the experimental data mea- 
sured at ENS by empty squares for an attractive Fermi gas at 
unitarity |4|,|7||, which agree extremely well with the prediction 
from quantum virial expansion. 



virial expansion up to the third order (VE3) and second 
order (VE2), respectively. For comparison, we also show 
the ideal gas result by the thin dot-dashed line. 

For a strongly attractively interacting Fermi gas, we 
observe that the quantum virial expansion is valid down 
to the degeneracy temperature Tp, where the predictions 
using the second-order or third-order expansion do not 
greatly differ. We note that our prediction of the third 
virial coefficient of a unitarity Fermi gas, Afog'^^ ^iom — 
—0.35501298, was experimentally confirmed to within 5% 
relative accuracy in the most recent thermodynamic mea- 
surement at ENS by Nascimbene and co-workers [7|- 

However, for a strongly repulsively interacting Fermi 
gas, the applicability of the quantum virial expansion 
is severely reduced: it seems to be applicable only for 
T > STp. Below this characteristic temperature, the dif- 
ference between the second-order and third-order predic- 
tion becomes very significant. This is partly due to the 
large absolute value of the third virial coefficient, sug- 
gesting that in this case the virial expansion converges 
very slowly. 



B. Harmonically trapped equation of states 

We finally present in Figs. 6 and 7 the high- 
temperature expansion prediction for the equation of 
state of a harmonically trapped Fermi gas in the strongly 
interacting regime. Due to the significantly reduced virial 
coefficients, the virial expansion in a trap has much 
broader applicability. For a strongly attractively inter- 
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Figure 7: (Color online) Entropy per particle S/{NEf) as 
a function of reduce temperature T/Tf for a trapped Fermi 
gas with infinitely attractive and repulsive interactions. The 
others are the same as in Fig. 6. 



acting Fermi gas, it is now quantitatively applicable down 
to O.STp, as confirmed by the precise experimental mea- 
surement at ENS (empty squares) [3, 0|. At the same 
time, the virial expansion for a strongly repulsively in- 
teracting gas seems to be qualitatively valid a,t T > Tp. 
At this temperature, the energy of the repulsively inter- 
acting Fermi gas is only marginally higher than the ideal, 
non-interacting energy. Considering the large energy dif- 
ference between a fully polarized Fermi gas and a non- 
polarized Fermi gas (i.e., at the order of NEp), we con- 
jecture that a strongly repulsively interacting Fermi gas 
does not have itinerant ferromagnetism in the tempera- 
ture regime where the quantum virial expansion theory 
is applicable. 



VI. CONCLUSIONS AND OUTLOOK 

In conclusion, we have presented a complete set of ex- 
act solutions for three attractively interacting fermions in 
a harmonic trap, with either positive or negative scatter- 
ing lengths. Firstly, we have outlined the details of our 
previous studies on the quantum virial expansion [6|, in 
particular the method for calculating the third virial co- 
efficient which was recently confirmed experimentally. In 
addition, we have opened up the previously unexplored 
repulsively interacting regime, and have presented a few- 
body perspective of itinerant ferromagnetism. We have 
also studied the high-temperature thermodynamics of a 
strongly repulsively interacting Fermi gas, by calculating 
its second and third virial coefficients in the unitarity 
limit. 

On the positive scattering length side of a Feshbach res- 
onance, a repulsively interacting Fermi gas is thought to 
occur by excluding all the many-body states which con- 
tain a molecule-like bound state for any two atoms with 
unlike spins. Strictly speaking, this is a conjecture which 



stems from a two-body picture. We have examined this 
conjecture using the exact three-fermion energy spectrum 
near the resonance. We have found some horizontal en- 
ergy levels that may be identified as the energy spectrum 
of three "repulsively" interacting fermions, as well as some 
vertical energy levels involving a tightly-bound molecule. 
However, many avoided crossings between horizontal and 
vertical levels make it difficult to unambiguously identify 
the energy spectrum of a repulsive Fermi system. 

For three "repulsively" interacting fermions in a har- 
monic trap, we have shown that close to the resonance, 
the ground state energy is higher than that of three fully 
polarized fermions. This is an indication of the existence 
of itinerant ferromagnetism in a trapped strongly repul- 
sively interacting Fermi gas. We have also considered 
the possibility of itinerant ferromagnetism at high tem- 
peratures. We have found that it does not exist in the 
regime where a quantum virial expansion is applicable. 
This gives an upper bound (^ Tp) for the critical ferro- 
magnetic transition temperature. 

Our high-temperature equations of state of a strongly 
repulsively interacting Fermi gas have a number of po- 
tential applications. We anticipate that these results 
can provide an unbiased benchmark for future quantum 
Monte Carlo simulations of strongly repulsively interact- 
ing Fermi gases at high temperatures [33 35] , using either 
hard-sphere interatomic potentials or resonance interac- 
tions. These results are also directly testable in future ex- 
perimental measurements, as inspired by the most recent 
thermodynamics measurement at ENS that have already 
confirmed our predicted second and third virial coeffi- 
cients for strongly attractively interacting fermions Q- 
Our exact three-fermion solutions in 3D harmonic traps 
will also be useful for understanding the dynamical prop- 
erties of strongly interacting Fermi gases at high temper- 
atures, by applying a similar quantum virial expansion 
for the dynamic structure factors [36,] and single-particle 
spectral functions [37|. 

These exact solutions of three interacting particles can 
be generalized to other dimensions, by adopting a suit- 
able Bethe-Peierls boundary condition for the contact in- 
teractions. Of particular interest is the case of two di- 
mensions, where the reduction of the spatial dimension- 
ality increases the role of fluctuations and therefore im- 
poses severe challenges for theoretical studies. The three- 
body solutions in 2D and the resulting high-temperature 
equations of state of strongly interacting systems will 
be given elsewhere, and provide a useful starting point 
to understanding more sophisticated collective phenom- 
ena such as the Berezinsky-Kosterlitz-Thouless transition 
and non-Fermi-liquid behavior. 

Note added: On finishing this manuscript, we are 
aware of a very recent work by Daily and Blume [38|, in 
which the energy spectrum of three and four fermions has 
been calculated using hyperspherical coordinates with a 
stochastic variational appoach. Our exact results are in 
excellent agreement with theirs when there is an overlap. 



13 



Acknowledgments 



where 



This work was supported in part by the ARC Centre 
of Excellence, ARC Discovery Project No. DP0984522 
and No. DP0984637, NSFC Grant No. 10774190, and 
NFRPC (Chinese 973) Grant No. 2006CB921404 and 
No. 2006CB921306. 



Appendix A: Calculation of Cn„' 

In this appendix, we outline the details of how to con- 
struct the matrix element C„„' in Eq. (|23p . which is 
given by, 

Cnn' ^ j P^dpRnl (p) Rn'l (f ) «'(^P; i'/.n'), (Al) 


where 
Rni [p) = 



2n! 



T{n + l + 3/2) 



pl^-p^/2^(l^ + l,2) ^p2^ ^ (A2) 



is the radial wave function of an isotropic 3D harmonic 
oscillator and the two-body relative wave function is 

< - n~^l,n')U{-Ui,,, ^, ^p2)^^p(_3^2)^ (A3) 

Here, for convenience we have set d = 1 as the unit of 
length. Ln is the generalized Laguerre polynomial 

and U is the second Kummer confluent hypergeometric 
function. A direct integration for C„„' is difficult, since 
the second Kummer function has a singularity at the ori- 
gin. The need to integrate for different values of vi_n' also 
causes additional complications. 

It turns out that a better strategy for the numerical 
calculations is to write, 



^\ 



rel 
26 



E 

fe=0 



1 



r(fc + 3/2) V3 ' 



(A4) 



by using the exact identity, 



r(-^)[/(-^,|,x2) = ^ 



oo j-1/2 



L^ {-') 



(A5) 



fe=0 



Therefore, we find that 
1 



Gji 



E- 

A:=0 



Vl, 



r(fc + 3/2)^; 



2fc! 



(A6) 



C: 



nji'k 



p^dpRni (p) R 




(A7) 



can be calculated to high accuracy with an appropriate 
integration algorithm. In checking convergence of the 
summation over fc, we find numerically that for a cut-off 
JT-max (i-e.; n, n' < rimax), Clin'k vanishes for a sufficient 

large k > fcmax '^ 4ninax- 

In practical calculations, we tabulate C^„//j for a given 
total relative angular momentum. The calculation of 
Cnn' for different values of i>in' then reduces to a simple 
summation over k, which is very efficient. Numerically, 
we have confirmed that the matrix C„„' is symmetric. 

I.e., Lynn' — L^n'ri- 



Appendix B: Calculation of Si,„ 

The calculation of s;.„ seems straightforward by using 
the Bethe-Peierls boundary condition in hyperspherical 
coordinates (pS)) . However, we find that numerical accu- 
racy is low for large n and I due to the difficulty of calcu- 
lating the hypergeometric function 2F1 accurately using 
IEEE standard precision arithmetic. We have therefore 
utihzed MATHEMATICA software that can perform an- 
alytical calculations with unlimited accuracy. For this 
purpose, we introduce As;.„ = s;_„ — si^n- After some 
algebra, we find the following boundary condition for 
t = As,,„/2, 



sin [irt) — \ — 



\n+l 



7:{-iy^'rin + l + l + t) 



3 2^r{i + l)r{n+i + t) 



fit), (Bl) 



where we have defined a function 



f{t)= 2F1 



t,n + l + l + t,l + 



3 1 

2' 4 



(B2) 



The above equation can be solved using the MATH- 
EMATICA routine "FindRoot", by seeking a solution 
around t — 0. It is also easy to write a short program 
to solve Eq. (jBip continuously for n < rimax = 512 and 
I < Imax = 512. In a typical current PC, this takes sev- 
eral days. The results can be tabulated and stored in a 
file for further use. 
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